Particle-in-cell simulations of particle energization from low Mach number fast mode 

shocks 
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Astrophysical shocks are often studied in the high Mach number Hmit but weakly compressive 
fast shocks can occur in magnetic reconnection outflows and are considered to be a site of parti- 
cle energization in solar flares. Here we study the microphysics of such perpendicular, low Mach 
number coUisionless shocks using two-dimensional particle-in-cell (PIC) simulations with a reduced 
ion/electron mass ratio and employ a moving wall boundary method for initial generation of the 
shock. This moving wall method allows for more control of the shock speed, smaller simulation box 
sizes, and longer simulation times than the commonly used fixed wall, reflection method of shock 
formation. Our results, which are independent of the shock formation method, reveal the preva- 
lence shock drift acceleration (SDA) of both electron and ions in a purely perpendicular shock with 
Alfven Mach number Ma ~ 6.8 and ratio of thermal to magnetic pressure /? = 8. We determine the 
respective minimum energies required for electrons and ions to incur SDA. We derive a theoretical 
electron distribution via SDA that compares to the simulation results. We also show that a modified 
two-stream instability due to the incoming and reflecting ions in the shock transition region acts as 
the mechanism to generate coUisionless plasma turbulence that sustains the shock. 



I. INTRODUCTION 

Solar flares convert magnetic energy into flow and par- 
ticle energy via magnetic reconnection (see e.g., Priest 
and Forbes(2002)[l|, Zharkova et aZ. (201 1)0, and ref- 
erences therein). The outflows from such reconnection 
sites can exceed the fast magneto-sonic speed. Unlike 
the inflows, the outflows from reconnection sites are flow 
dominated and the ratio of thermal to magnetic pressure 
(= /3) exceeds 1. Analytic predictions Q and numerical 
simulations in which an obstacle is present [J, M, reveal 
the presence of low Mach number fast shocks in these re- 
connection outflows. In the standard geometry of a solar 
flare, such "termination" shocks form, where the down- 
ward directed outflow interacts with the magnetic loop 
formed from previously reconnected field lines. 

CoUisionless termination shocks have been invoked in 
a number of phenonienological solar flare models and 
may be associated with specific observational features. 
Mann et a/. (2006) (2009) [i Q suggested shock drift ac- 
celeration(SDA) as a mechanism of energetic electrons 
up to ~ MeV duringsolar flares. Hard X-ray emission 
at loop top locations f8| may also be associated with such 
shocks. Decker and Vlahos(1986) 9] carried out test par- 
ticle simulations on the SDA in solar flare shocks. Guo 
and Giacalone(2010)^J studied a solar flare shock using 
a hybrid simulation in the presence of turbulent mag- 
netic fields, where electrons are efficiently accelerated 
by multiple refiection. The hybrid simulations treat the 
ions as particles and the electrons as a fluid. However, 
to our knowledge, fully-kinetic simulations, where both 
shock fornration and particle acceleration are modeled 
from first principles in the regime of low-Mach-number 
and 13 > 1 have not been reported. 



Particle-in-Cell (PIC) simulations have been used to 
study high-Mach-number coUisionless shocks and parti- 
cle acceleration in ion-electron plasmas |10l - ll5| but Low- 
Mach-number shocks are less widely studied. Gargate 
and Spitkovsky(2012)[l7] investigated a low Mach num- 
ber shock as a subset of cases in a parameter survey using 
a hybrid code, where both diffusive shock acceleration 
(DSA) and SDA of the ion acceleration were observed, 
with the latter becoming more important as the shock 
becomes more perpendicular. 

Here we report results from two-dimensional (2D) full- 
PIC simulations of perpendicular shocks in the regime of 
low-Mach-number (< 3) and high plasma (3{> 1). The 
motivation is to study the microphysics of shock forma- 
tion and particle acceleration in the perpendicular shocks 
relevant to solar flares. Perpendicular shocks are chosen 
for their relevance to the shocks that emerge in the 2-D 
reconnection outflows 5] and for the first-step toward the 
study of more general quasi-perpendicular shocks. As we 
describe later, the simulations reveal both electron and 
ion acceleration via SDA. We also find that the modified 
two-stream instability by the interaction of incoming and 
reflecting ions in the shock transition region is the likely 
turbulent dissipation mechanism that sustains the coUi- 
sionless shock. 

Most previous PIC simulations of coUisionless shocks 
use a reflecting wall boundary condition where plasma 
flow reflects off a rigid wall to form a shock. In that 
case, the simulation frame is fixed to the downstream 
rest frame. The shock moves away from the reflecting 
wall at a speed Vg, the downstream flow velocity in the 
shock rest frame. The simulation time is then limited to 
Lx/Vs where Lx is the simulation box size in the direction 
of the shock propagation. 



In contrast our simulations use a moving wall bound- 
ary condition, first introduced in Langdon et al.^^. It al- 
lows control of the downstream flow velocity and thus the 
shock speed in the simulation frame. By slowing down 
the shock speed in the simulation box, smaller boxes can 
be used for the same simulation time. We have imple- 
mented this boundary condition in 2D and find all prop- 
erties of the generated shocks are similar to those with 
the reflecting boundary condition when differences be- 
tween the reference frames are accounted for. 

The rest of the paper is organized as follows. The 
simulation setup, including the moving wall boundary 
condition, is described in section |lTl The shock proper- 
ties and particle acceleration are described in section Hill 
Discussion and summary are given in section ITVl 



II. SIMULATION SETUP 
A. Basic setup 

We use the relativistic fuh PIC code OSIRIS pj,] to 
study the formation of low Mach number fast perpen- 
dicular magnetosonic shocks and the particle accelera- 
tion therein. Following Refs.j5|] and ^2011, typical param- 
eters of solar flare reconnection outflows are chosen as 
the upstream conditions for our shock. (Hereafter, "up- 
stream" is defined respect to the fast shocks we study 
herein, not upstream of a reconnection site.) In par- 
ticular, we use a plasma density n = 5 x 10^ cm^'^, 
electron and ion temperatures T^ — Ti = O.SkeV, and 
the perpendicular magnetic field strength _B = 5G with 
P = 8TTn{Te + T,)/B^ = 8.05. A reduced ion/electron 
mass ratio of mi/rrie — 30 is used to reduce the required 
computational resources. The Alfven Mach number is 
chosen to be Ma = Vi^/Airrrnn/B — 6.79, which im- 
plies an upstream plasma flow velocity in the shock rest 
frame, Vi — 0.0274c, where c is the speed of light, for 
mi/rrie = 30. For a real proton-electron plasma, Vi would 
be reduced by a factor of ^30/1836 = 0.128 for the same 
Ma. The super- fas t-magnetos onic Mach number M sat- 
isfies M = Ma/\JI + (5/6)/3 = 2.45. The ratio of the 
electron cyclotron frequency to the electron plasma fre- 
quency is Q,ce/0Jpe = 0.02207. 

With these upstream values of Ma and /3, the Rankine- 
Hugoniot relation[21j for perpendicular shocks gives the 
compression ratio r = 2.15, 
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where the lower indices 1 and 2 represent the upstream 
and the downstream, respectively, in the shock rest 
frame, and we used the adiabatic index 7 = 2 for 2D. 

A rectangular simulation domain in the xy plane is 
used. A uniform external magnetic field oi B = 5G is 
set out of the simulation plane (along the z-axis) and 
a uniform external i?-field(= VdB/c) is set up along 



the y-axis. (An alternative in-plane configuration with 
By and E^ has also been used to help identify insta- 
bilities responsible for dissipation and has yielded sim- 
ilar results.) The simulation box is initialized with a 
Maxwellian ion-electron plasma drifting to the right with 
Vd = 0.021c and T^ = Tt = 0.5keV, where Vd is set to 
a smaller value than the upstream speed Vi{= 0.0274c) 
in the shock rest frame considering the shock speed trav- 
eling to the left direction. A new plasma of the same 
distribution is constantly injected from the left bound- 
ary [x — 0) throughout the simulation. The simulation 
box sizes are L^ x Ly ~ 340c/a;pe x 40c/wpe. These 
correspond to 62c/u!pi x 7.3c/u!pi, where tUpi is the ion 
plasma frequency, for mi/rrif, — 30. The grid size used 
is dx = dy = 0.08c/wpe and the time step used is 
dt — O.OSG/wpe. For diagnostic purposes, a small pop- 
ulation of electrons and ions spatially localized within a 
circular region, is set as a separate species for which the 
particle information is stored more frequently to track 
detailed trajectories over time. For each particle species, 
25 particles per cell are used. A linear current deposition 
scheme is used for all simulations in this paper. 

A periodic boundary condition is used in the y- 
direction for both particles and fields. For fields, an open 
boundary condition is used in the x-direction. Particles 
that reach a; = are re-injected into the box with the 
initial drifting Maxwellian distribution. Ai x ~ Lx, & 
moving wall boundary condition is adopted, as described 
in the next subsection. 



B. The moving w^all boundary condition 

In general, a refiecting wall moving in the direction of 
the flow[lg| can force the plasma flow velocity at the wall 
to be an arbitrary predetermined velocity by selectively 
reflecting particles with certain velocities. Here we im- 
plement a moving wall boundary condition at a; = L^^ to 
force the flow velocity there to be close to the downstream 
velocity measured in the shock frame. The wall is essen- 
tially treated as an infinitely massive slab moving with 
velocity Wwaii and particles that catch up to it rebound 
specularly off the wall in the wall rest frame. Proper rel- 
ativistic momentum transformations are applied to ob- 
tain the particle velocity in the simulation frame after 
rebounding. Between times t„ and i„ -I- At, the wall will 
have "moved" a distance KvaiiAi from the right domain 
boundary L^- Particles will also be moving beyond the 
domain boundary during this time. Particles that can 
reach the wall and rebound quickly enough to return to 
the simulation domain during At are kept in the box. 
Those which do not are removed from the system. At 
the right boundary, this procedure forces the bulk flow 
velocity to be Kvaii- 

One immediate consequence of this implementation is 
that it cannot be used to probe shocks with compression 
ratios of r = 2 or less. This can be seen by noting that a 
necessary condition for the plasma to return to the simu- 



lation domain is that its updated velocity must be nega- 
tive in sign. In one dimension, a non-relativistic particle 
with a velocity Vp will rebound with a velocity 2V„aii — Wp. 
This must be negative for the particle to remain in the 
box. Initially, Vp « Vd where Vd is the plasma flow veloc- 
ity from the left boundary. In the shock rest frame, Vd 
is the upstream velocity and Vwaii should be the down- 
stream velocity, Vivaii = Vd/r. The condition for rebound 
back into the box then becomes 2Vd/r — Vd < 0, where 
r = n2/ni is the compression ratio across the shock. For 
r < 2 the rebound condition cannot be satisfied for purely 
elastic collisions with the wall. 

In practice we do not set the wall velocity to be ex- 
actly the downstream velocity in the shock frame but 
set it so that the shock is slowly propagating back into 
the upstream. This is necessary to allow a large enough 
downstream region to be generated wherein particles may 
travel several ion gyro-radii to undergo acceleration. We 
find that we can control the shock velocity in the sim- 
ulation frame by changing Vwaii and have tested that 
the properties of these shocks are essentially the same as 
those generated from a stationary reflecting wall bound- 
ary. 
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III. RESULTS AND ANALYSIS 

A. shock structure 

In Figdl we plot the ratios of the density, velocity, 
and magnetic field [n(a;)/ni, Vi/Vx{x), and Bz{x)/Bi], 
the momentum phase space distributions oipxX and PyX, 
the y-averaged flow velocity profiles Vx (x) and Vy (x) , and 
the temperature profiles Tx{x) and Ty{x), for the ions 
(the left column) and the electrons (the right column) 
at i = 28476/ajpe. Here, the momenta are defined as 
Px — Tvx and Py — Tvy, where F = I/-\/l — (w/c)^. The 
ratio curves [(a) and (b)] are calculated in the shock-rest 
frame using the Lorentz transformation and the other 
plots are computed in the simulation frame. 

In FiglU the shock front is at a: « 150c/wpe and moves 
to the left with a shock speed in the simulation frame of 
6.4 X 10~^c. An oscillatory pattern in the downstream 
properties can be observed, which indicates weak dissi- 
pation in this low-Mach-number shock. The Rankine- 
Hugoniot condition in Eq.([T|) gives r = 2.15 for our up- 
stream parameters. In the simulation, the compression 
ratio is r ~ 2.8 near the shock front and relaxes to r ^ 2.1 
in the far downstream with a weakly damped oscillatory 
pattern (FigH^ and b). In the downstream, the elec- 
trons are thermalized isotropically (Fig[TJi, f and j) but 
the ions are heated slightly more in the y-direction near 
the shock front (Fig[TJ; and i). The electrons are heated 
to Ta « 1.4keV and T2/T1 « 2.8 (Fig[T]), which indicates 
that the electrons are mainly heated by adiabatic com- 
pression, {T2/Ti)adia = in2/niy~'^ = r^'^^ = r, where 
the adiabatic index 7 is 2 in 2D. The electron tempera- 
ture in the downstream slightly increases as the compres- 



FIG. 1; (Color online) The ratios, n/ni, Vi/V^, and B^/Bi, 
momentum distribution of p^, and py, y-averaged flow velocity 
of Vx and Vy, and temperature Tx and Ty, for ions(left col- 
umn) and electrons (right column) at i = 28476/a;pe. (a) ~ (6) 
are calculated in the shock-rest frame and the other plots are 
obtained in the simulation frame. 



sion ratio relaxes to ~ 2.1 (Fig[Tt). The ions are heated 
to a value much higher than the electrons near the shock 
front (namely T2/T1 « 3 ^ 5 in FigHJ), highlighting the 
more substantial role of non-adiabatic particle energiza- 
tion for the ions than the electrons. 

The flow velocity in the downstream is mainly in the 
x— direction, osciUating around Kvaii(= 0.0052c) (Fig[T^ 
and h). The flow velocity in the i/-direction oscillates 
around zero in the downstream (FiglT^ and h) due to 
the Ex X B-drift. The damped oscillatory pattern in the 
downstream results from weak turbulent dissipation in 
the shock transition region via the modified two-stream 
instability ^,,23;]. This instability arises from the incom- 
ing and reflecting ions at the shock front. The ions are 
reflected in a thin region, 140 < x < 155{c/ujpe) (Figdt), 
due to a potential jump, eA$, across the shock front 
given by the electron momentum cquation|24l |25|| , 
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where the ion drift Viy is omitted. 

Figure [5] shows the y-averaged Ey, Ex, and the electric 
potential energy e^{x) at t — 28476 /ujpe in the shock 
rest frame. Ey is positive and approximately constant 
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FIG. 2: j/-averaged Sj,, -Bj,, and the potential energy e4>{x) 
at t = 28476/a;pe in the shock rest frame. In (c), the shock 
transition region is indicated by the dotted lines. 
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FIG. 3: (Color online) (a)(c)ion and (d)(d)electron distribu- 
tions in the shock transition region, 140 < a; < 155(c/a;pe) in 
the simulation frame. In (c) and (d), the distributions fit into 
Maxwellian distributions (dashed lines). 



across the shock (Fig[2^), but a negative E^ causes a 
potential barrier for the incoming ions at the shock front 
(Figl^b and c). The potential energy is e^{x) « 3kcV [a 
bit larger than the 2.3keV from Eq.([2])] and can reflect 
the low energy tail of the incoming ions, which have an 
average drift energy in the shock rest frame of 5.74keV 
{Vi = 0.0274c). 



1. Modified Two-Stream Instability as the Source of 
Dissipation for Shock Sustenance 

Figure [3] shows the electron and ion distributions in 
the shock transition region, 140 < a; < 155(c/wpe). The 
ions have a bump-on-tail distribution, with 22% of the 
ions reflected (Figl3^ and c). A high energy component 
can also be observed moving in the positive-?/ direction 
(FiglSJi.), resulting from SDA, which will be discussed 
later. The electrons are essentially isotropic and drift 
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FIG. 4: (Color online) (a)~(d): Numerical solutions of 
the modified two-stream instability under B = 6.8G, Ti = 
0.75keV, and Te = 0.85keV (solid lines), and different param- 
eters (dotted and dashed lines), (e) and (f): Fourier spectra 
of Ex and Ey fields from the simulation at f = 28476/a;pe . 



with Vex = 0.0085c(Figl3)D and d). We use the distribu- 
tion in Fig|3l[c) and (d) for a linear stability analysis of 
the modified two-stream instability(MTSI) at the shock 
transition region in the electron drift rest frame to as- 
sess whether the associated dissipation is consistent with 
what is needed to sustain the shock. 

For the stability analysis, we also use the following ini- 
tial parameters extracted directly from the simulation: 
The magnetic field is B^ = 1.36.Bi = 6.8G. The electrons 
are magnetized and have a Maxwellian distribution with 
a temperature of To = 0.85keV. The ions are assumed 
to be non-magnetized and have drifting Maxwellian dis- 
tributions. In the electron rest frame, the drift veloc- 
ities and densities for the incoming and reflecting ions 
are 14in = 0.0075c, V^re = -0.0165c, njn = 0.78ni, 
Uro = 0.22ni. Both incoming and reflecting ions have 
a temperature of Txin = T^-cc = 0.75keV. 

The dispersion relation for the MTSI with k = kx is|26| 
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where Vsthi= \/Ts/ms) is the thermal velocity of species 
s{= e,in,re), Xe = k'^v^th/^le^ ^m(Ae) is a modified 
Bessel function of the 2nd kind, £_s — {i^ — kVxs)/V2kvsth, 
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FIG. 5: (Color online) field fluctuation, 5Ex, 5Ey, and SB^ at 
i = 28476 /ojpe. In (c), SEx/SExo is plotted and the red line is 
exp[7(a; — xo)/Vi], where 7 = 4 x 10~'^LUpe from the modified 
two-stream instability, xo = lAQc/tOpe, and Vi = 0.0274c. 



and Z(^) is the plasma dispersion function. We then nu- 
merically solve Eq.Q, using a fractional polynomial ap- 
proximation of the Z{^) function[27], the Zenkins and 
Traub algorithm ^28i, 2M for complex polynomial root 
finding, and the Muller method[30] to obtain accurate 
numerical solutions. 

Figure|4]shows numerical solutions of Eq.Q for the pa- 
rameters given in the paragraphs above (solid lines) and 
also for slightly different parameters (dotted and dashed 
lines). The maximum growth rate in the solid line is 
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at fc = 0.2ujpe/c (Figgli). The growth 



rate increases as the magnetic field increases (Fig|3^) 
and decreases as the ion and/or electron temperatures 
increase (FigHJ: and d), but is less sensitive to Tg than 
to Ti. In FigllKb), the real frequency is negative, giving 
a phase velocity comparable to the drift velocity of the 
reflected ions in the electron rest frame. We compare 
the numerical solutions with the modes observed in the 
simulation: In Figllje) and (f), the Fourier spectra of 
Ex and Ey fields from the simulation are plotted. An 
electrostatic mode is observed at k^ = 0.2ujpe/c in E^, in 
agreement with the MTSI dispersion relation. 

Figure[5]shows the fluctuating fields along the x-axis at 
t = 28476/ u!pe for SE^, SEy, and SBz, where we define 
SA = ^{{A — (A))^) and (..) is the ensemble average. 
Here we use the y-average as the ensemble average. In 
(a) and (b), SE^ and SEy increase at the same rate in 
the shock transition region. In (c), we plot the ratio 
of SEx to 6Exo, where SE^o is the fluctuating field in 
the upstream. In (d), the magnetic field fluctuation is 
negligible compared to the electric field fluctuation since 
the MTSI has electrostatic modes. The red line in (c) 
is the evolution curve of the fastest growing mode in the 
MTSI, exp[7(a; — xo)/Vi], where the growth rate 7 — 
4. X lO-'^Wpe and the shock speed Vi = 0.0274c. The 
MTSI begins to operate at a: = xq = UQc/ojpe and ends 
before x = 160c/iOpe during the shock transit time. At ^ 
Aa;/Vi ~ 730/ujpe, and is enough to excite the electric 
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FIG. 6: (Color online) (a) specific entropy s{x) for ions (b) 
F{x) for ions (c) As = s(x) - si (red), AF/mV^i = {F{x) - 
i^i)/7iiKi (blue). As - AF/niVxi (black) for ions. 



field fluctuation observed in the shock transition region. 
Whether this level of SE is sufficient to generate the 
entropy increase required for the shock is an interesting 
question. The entropy equation can be written as by 
acting / dv{l + ln(/)) on the Vlasov equation |21||31| 



d{ns)/dt + d{nVxs)/dx = d/dx / dv'<(/)ln(/) 
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where s{x) is the specific entropy (entropy/particle), 
six) = -J d^r{{f{x, v)) ln(/(a;, v))}/ / d^r{f{x, v)), n = 
/ dv(/), T4 = (1/n) / d-vvxif), and w^ = -y^, - T4. In the 
shock rest frame {d/dt = 0), Eq.(|3]) can be written as 
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Here F{x) ^ j dVv'x{f)hi{f) as seen in Ref.|31|. 

In principle, one could measure every term in Eq.([5]) 
and show definitively whether this level of 5E is sufficient 
to generate the entropy increase required. However, the 
term on the right-hand side, the entropy change due to 
fluctuating fields, involves 5f and is difficult to evaluate 
from the simulation. Here, we measure the terms on the 
left-hand side of Eq.([5)). Figure [6^ a) shows the specific 
entropy for the ions measured from the raw particle data 
in the simulation. The specific entropy jump between 
the upstream and downstream is ^ 0.4. The Sackur- 
Tetrode equation for the entropy of an ideal gas gives 
a similar result for the entropy jump in this simulation. 
As = {S2-Si)/N = ln{(l/r)(r2/ri)'*/2| ^ q 35^ ^t^^^^ 

d is the degree of freedom (d = 2 for 2D), r = 2.1, and 
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FIG. 7: Simulation result of the energy distribution, /(e) = 
{l/Ntot)dN{e)/de, in the downstream rest frame of 145 < 
X < 340{c/ujpe) at i = 28476 /ojpe. (a) Ion energy distribution 
and fittings into Maxwellian temperatures with Ti = 1.5(dot- 
dashed), 2.6(dotted), and 1. 6 (keV) (dotted lines), (b) Elec- 
tron energy distribution and fitting into Maxwellian temper- 
atures with Ti = 1.4keV(dashed line) in < e < 10 and 
T = 1.9keV(dotted line) in e > 13(keV). 



T2/T1 — 1.5/0.5. Figure [6ljb) shows F{x) for the ions in 
Eq.(l5]). When (/(«)) is an even function about Vx such 
as a drift-Maxwellian distribution, F{x) vanishes. We 
see non-zero F(x) in the transition region. In FiglHl^c), 
we plot As = s(x) — si (red line), AF/niVxi — {F{x) — 
Fi}/niVxi (blue line), and As — AF/niVxi (black hne) 
for the ion. Therefore, we conclude that the fluctuating 
fields excited by the MTSI are necessary for the entropy 
creation throughout the downstream region. Whether 
they are sufficient remains an open question. 

Another possible instability responsible for shock for- 
mation comes from diamagnetic currents on the shock 
front [2 Ij. and would have a mode with ky. This mode, 
however, is not observed in either E^- or i?y-spectrum 
as seen in Fig|3t^e,f). This is supported by another per- 
pendicular shock simulation where the magnetic field is 
initiated in the simulation plane and the fcy-mode is pre- 
cluded but the same shock structure is seen. 



B. Particle heating via shock drift acceleration 

1. Dynamics of ions and electrons incurring SDA 

Shock surfing acceleration (SSA)[l3, [H [H, HI, HI, 
whereby particles are trapped in the solitary wave struc- 
ture excited by the Bunemann type instability at the 
shock front and then accelerated by the convective Ey 
field along the shock, is inefficient for low Mach perpen- 
dicular shocks. We find no evidence of SSA in our simu- 
lation. However, when particles pass through the shock 
transition region, the magnetic field jump at the shock 
front allows particles to experience shock drift accelera- 
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FIG. 8: A typical ion tracking experiencing SDA (in the left 
column) and not experiencing SDA (in the right column). In 
(a), we plot the positions of the shock front, the 2nd, 3rd, 
and 4th peaks traveling to the —x direction in the simulation 
frame. 



tion (SDA)[i,[3i-l43. SDA results from the fact that the 
Ey is constant across the shock but magnetic field jump 
gives different gyro-radii of a particle at two sides of the 
shock. This gives rise to a net drift along the y-axis, Jy, 
and a net gain energy of 5e = eEySy per particle. 

In our simulation, both ion and electron heating via 
SDA are observed. Figure [7] shows the normalized ion 
and electron distributions in the downstream rest frame, 
/(e) = {l/Ntot)dN{e)/de, in the downstream region of 
145 < X < 340(c/a;pe) at t = 28A76/ujpe. Both distribu- 
tions have a low energy regime that corresponds to adia- 
batic heating(e < 5keV for the ion and e < lOkeV for the 
electron) and a high energy regime due to SDA. The ion 
distribution fits into a multi-temperature Maxwellian dis- 
tribution with temperatures of T^i = 1.5keV in < e < 
5(keV), r,2 == 2.6keV in 5 < e < 17(keV), and T,3 = 
1.6keV in e > 17keV (FigEJi). The electron distribu- 
tion fits into a two-temperature Maxwellian distribution 
with temperatures of Tei — 1.4keV in < e < lO(keV), 
Te2 = 1.9keV in e > 13(keV) (FigEb). The minimum en- 
ergy where SDA is effective has different origins for the 
ions and electrons. The difference mainly comes from the 
fact that the magnetic moment is conserved for the elec- 
trons but not for the ions in the shock transition region. 

We discuss SDA for the ions first. Figure [5] shows typ- 
ical tracks for an ion experiencing SDA (the left column) 
and for one not experiencing SDA (the right column) 
from the simulation. Figure E^a) shows the particle's 
x-coordinate vs time, along with positions of the shock 
front and the subsequent 3 compression peaks identifiable 
in Fig[T]Ja). When the particle meets the shock front, it 
turns back toward the upstream. The trajectory in the 
xy-plane in Fig[5]Jb) shows the particle drifting up along 
the y-axis with a larger gyro-radius after turning back. 
Accelerated by the Ey field, the kinetic energy of the 
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FIG. 9: (a) A schematic view of the ion trajectory in the 
flow rest frame in the shock transition region where the shock 
front travels to the left with speed V!,. (b) The shaded re- 
gion represents the ions experiencing SDA in the flow rest 
frame. The minimum velocity to experience SDA is indicated 
by Smin(= v/V!,). 



ion increases from 4keV up to 16kcV after encountering 
the shock front [FigHfc]. This particle later re-crosses the 
shock front and passes through the secondary compres- 
sion peaks where it loses some energy, reaching the right 
boundary with an energy of 8keV. The net energy gain of 
the particle results from SDA at the encounter with the 
shock front, despite some energy loss during interaction 
with the secondary peaks. In Figl8l[d)-(f), a different ion 
entering the shock front with a different angle and energy 
does not meet the required conditions (discussed below) 
to turn back to cross the shock front again and the ki- 
netic energy of the ion decreases to IkeV after reaching 
the right boundary. 

Whether or not an ion gains energy via SDA in a per- 
pendicular shock depends on its incident speed and an- 
gle of incidence at the shock front, and this is similar to 
the case of the electrons in a relativistic perpendicular 
shock[33, |33- The minimum energy above which SDA 
is operative is best computed in a reference frame where 
the ion drift velocity nearly vanishes and the ion executes 
approximately pure gyro-motion. Because the upstream 
and downstream drift velocities are different, the refer- 
ence frame is chosen such that the shock front moves to 
the left with speed V^ = l/2(I/i + V2) = V^i(l + r)/{2r) 
(Figin^i'). Since the gyro-radius for an ion is larger in 
the upstream than in the downstream, an ion will only 
gain energy if it returns to the upstream side drifting up- 
ward. A necessary condition for the ion to return to the 
upstream side, as shown in FigEJa), is 



p(l - sine) > I/'Ai, wcos6i > -y;. 



(6) 



where p = v/Qc is the gyro-radius (w > is the par- 
ticle speed and Vie is the gyro- frequency) , and At = 
(37r/2 + 9)/Hc is the time between the shock front cross- 
ings. Equation (jG)) can be re-written as 
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FIG. 10: A typical electron tracking experiencing SDA (in the 
left column) and not experiencing SDA (in the right column) 
in the simulation frame. 



with a dimensionless variable v = v/Vg. In Figl^l^b), 
Eq.([7l) is plotted in the VxVy-'pla.ne where the ions in 
the shaded region satisfy Eq.(I7]). The minimum ve- 
locity for the ions that can gain energy through SDA 
is indicated by the shortest distance from the origin 
to the shaded area in FiglQlJb) and is measured to be 
"Smin ~ 1.38. The corresponding minimum energy is 
emin = l/2M(z;niiny/)^=5.04 keV, using the measured 
value of r = 2.8. This agrees reasonably well with the 
transition at emin ~ 5keV in Fig[71[a), measured in the 
downstream rest frame. 

We now describe electron SDA. In FigHUl we plot typ- 
ical tracks of an electron experiencing SDA (the left col- 
umn) and one not experiencing SDA (the right column) 
from the simulation. Unlike the ions in FigUl the elec- 
trons drift through the shock fFigflUk and d) without 
turning back. Since the electron gyro radius is small 
compared to the shock width, the electron magnetic mo- 
ment can be treated as a constant in the shock transition 
region. In addition to the drift in the x-direction, elec- 
trons also drift in the y-direction due to Ej; x B-drift and 
VB-drift. 

In Fig fTUT b). the electron drifts downward {~y axis) 
due to the VB-drift in the shock transition region around 
X — 260c/a;pe and its kinetic energy increases up to 20keV 
after encountering the shock front [FigflUb]. After leav- 
ing the shock front, the energy decreases to IkeV in the 
dB/dx < region and then increases again to 17keV once 
it encounters the 4th compression peak. In contrast, the 
electron in FigfTUTe') drifts upward {y axis) due to the 
Eix X B-drift (i?2: < in the shock transition region) and 
its energy is fluctuating between < e < 3 keV. 

For an electron to experience SDA, it must drift down- 
ward, which is possible when the VB-drift is larger than 
the Ea; X B-drift in the shock transition region where 
dB/dx > and -E^: < 0. This condition can be written 
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> 



2eB2 Ax B^ ' 



(8) 



where Ax is the width of the shock transition region. The 
minimum energy for the electrons to experience SDA is 
then 



Ax e\EJl:^x fr+l 



AB. 



r- 1 



, (9) 



where B^ and AB^ are taken as {r+l)Bi/2 and (r— l)i?i, 
respectively. With r = 2.8 and E^ — —0.0003 mcwpe/e 
for our shock, Cmin is estimated to be 3.24keV. The neg- 
ative Ex field in the shock transition region requires the 
electrons to have higher threshold energy for SDA. 



2. Electron Spectrum from SDA 

Now we discuss the electron spectrum due to SDA. 
Neglecting adiabatic heating, the electron energy change 
during one gyro-cycle in the shock transition region is 
i^esDA = sEySy, where Sy is the net drift distance from 
the VB-drift and the E^; x B-drift during one gyro-cycle 
and is given by 



Sy 



ec AB^ c\Ex\\ 27r 



zBl Ax B, n 



(10) 



Here, e(= mv^ 12) is the electron energy of gyro- motion. 
The rate of energy change from SDA is given by 



de ae — rj 
di ~ T ' 



(11) 



where r — 27r/i7c is the period of one gyro-cycle, and a 
and 77 are defined respectively as 

2TrmcViBiAB, 2TimcViBi\Ex\ ,_,^, 

a= — ^,»y= TT. • (12) 



eBl Ax ' 
The solution to Eq.dlTl) is 



(ei 



jexp 



at 

T 



Bl 



+ emin, (ei > Cmin), (13) 



where t\ is the particle energy before entering the shock 
transition region and emin is written as emin = ??/a using 
Eq.® and ^. 

If we assume that the particle enters the shock at a 
time t = with energy t\ then the probability for a 
single electron to leave the shock transition region at a 
time t is given by a Dirac delta function. 



P(t)dt = 8(t - T)dt, 



(14) 



where T is a characteristic time for escape and is given 
by T = Ax/ (14;), Ax is the width of the shock transi- 
tion region, and (Vx) is the averaged drift velocity. The 



associated probability of an electron to change its energy 
from ei to e is 

P(e, ei)de ^ S [e - {ei - e,nin)e"^/" - e,,i„] de. (15) 

The energy distribution in the downstream rest frame is 
obtained by 



/•oo 

/2(e2)=/ dei/i(ei)P(e2,ei), 



(16) 



where /i(ei) is the energy distribution in the upstream 
rest frame and given by {l/Ti)e~'^'-^'^'- . Then the energy 
distribution after SDA is a translated Maxwellian distri- 
bution with a temperature T2 = Tie"^'"^ for £2 > Cmin, 



/2(e2) 



y^gaT/r 



exp 



Ti 



r^f2 + (e"^/^-l).„in} 

(17) 



where aT/r is given by 

aT _ ViBi AB^ _ 8r(r - 1) 
~r " i32 (Vx) ^ {r + 1)3 ■ 



(18) 



For r = 2.8 in our shock, e"^/^ w 2.09 and T2 = 
y^gaT/r ^ i.04keV. 

If we include the adiabatic heating in Eq. (fTTj) , then the 
energy equation in the shock transition region becomes 



de _ a{e- emjn) 
dt ~ T 



C + i' 



(0 < t < T, ei > £.„;„), (19) 



where ^ = T/(r — 1) and the solution is given by 



e = e"*/^ f 1 + - 



.i+...„^e"«/^{i?J-^ 



-E, 



(-7(^ + ^))}]' (0 < ^ < T, 61 > e„,i„), (20) 



where Ei{s) = j'_^{e'i / q)dq. 

When the escape probability for an electron in the 
shock transition region is given by Eg.dTi)). we calculate 
the energy distribution in the downstream using Eq. (jl6p 
and (PO]) . We get a translated Maxwellian distribution 
with a temperature T2 = Tire"^'"^, 



/2(e2) 
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exp 



y^^gQT/r 






X E,, - 



-E,, - 



a^ 



, (£2 > £*), 

(21) 



where e, — e(t = T, ei = emin) in Eq. ([20| . aT/r is given 
by Eq.([TH]), and a^/r = 8r/(r -t- 1)^. For r = 2.8 in our 
shock, e, = 12.9keV and T2 = Tire"^/^ = 2.9keV. The 
energy distribution in the downstream for 62 < Cmin is 
given by pure adiabatic heating. 



72(62) = ^^rcxp 
rTi 



£2 

rTi 



, (£2 < emin). (22) 
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FIG. 11: Theoretical electron energy distribution /2(e) after 
SDA-f adiabatic lieating(solid line) and pure adiabatic heat- 
ing(daslied line) with r — 2.8 in (a). In (b), r — 2.1 for 
SDA-fadiabatic heating(solid line) and r — 2.8 for pure adi- 
abatic heating(dashed line) are taken. The simulation re- 
sult(thick solid line) (Fig (TJd) is plotted for comparison. 



In FigHH we plot Eq.dH]) (solid line), Eq. (gH) (dashed 
line) and the simulation result(thick solid line)[Figl7}D]. 
The temperature T2 = 2.9keV for e > e,(= 12.9keV) 
from Ea. (l2T]) . however, is larger than the T2 — 1.9keV 
from the simulation(Figl7)D). The difference between the 
theory and the simulation results from the oscillatory 
shock structure in the simulation where the high en- 
ergy electrons gain and lose energy via SDA by drifting 
downward(— y axis) in the region of dB/dx > and by 
drifting upward(-|-y axis) in dB/dx < 0, respectively. 
When an effective compression ratio, r — 2.1, is ap- 
plied to the electrons experiencing SDA, the correspond- 
ing temperature becomes Ti = T\re°''^ 1'^ = 1.9 keV for 
e > eH.(= 8.6keV), which is consistent with the simula- 
tion result. In FigfTTTb). we plot Eg. (|2T|) (solid line) with 
r = 2.1, Eg. ([^ (dashed line) with r = 2.8, and the sim- 
ulation result(thick solid line). 



C. Effects on spectra w^hen a realistic mass ratio of 
protons and elections is used 

Since our simulations have employed a reduce mass 
ratio of rrii/me = 30, it is natural to wonder what the 
energy spectra would be for a real mass ratio mijmf^ — 
1836 for fixed Mach number and plasma /3. The energy 
spectrum depends on the shock structure, in particular 
the changes of B^ and £'x across the shock. 

The change in B^ across the shock from B\ to rB\ 
and the shock electric potential jump eA$ in Eq.([5]) does 
not vary with the ion mass for fixed Mach number and 
plasma /3. Thus, the same fraction of ions as computed 
for our reduced mass ratio will reflect at the shock front 
for the realistic mass ratio case. The MTSI still operates 
and the growth rate is reduced by ^30/1836 for fixed 



M ach num ber and plasma /3. This implies that we need 
a/1836/30 times longer time to generate the shock in the 
electron plasma frequency (1/wpe) timescale. In addition, 
the electron energy distribution via SDA in Eg. (PT|) as 
well as the minimum energy for the electrons to incur 
SDA in Eq.(|9l) does not vary with the ion mass. 

In short, we expect the energy spectrum at least for 
electrons observed for our simulations with a reduced 
mass to be unchanged when a realistic mass is used. 



IV. CONCLUSIONS 



We simulated a purely perpendicular, low Mach num- 
ber fast mode shock relevant to the termination shocks 
in solar flare magnetic reconnection outflows. We em- 
ployed a 2D particle-in-cell code with a moving wall 
boundary condition for a reduced ion/electron mass ra- 
tio mi/rrie = 30. The moving wall method generates a 
slowly propagating shock compared to the more standard 
fixed reflection boundary method, and allows for smaller 
box sizes and more efficient use of simulation time. 

Both electrons and ions experience the shock drift ac- 
celeration(SDA) in our simulations. The transition en- 
ergy point from pure adiabatic heating to SDA measured 
in the downstream rest frame is given by ^ 5 keV for the 
ions and ~ 10 keV for the electrons. These values are well 
modeled by our theoretical analysis of SDA. The nega- 
tive Ex field in the shock transition region requires the 
electrons to have high threshold energy for SDA. The 
ion energy distribution in the downstream shows a tri- 
Maxwellian distribution and the electron energy distri- 
bution in the downstream shows a bi-Maxwellian distri- 
bution. We theoretically modeled the electron energy 
distribution via SDA with/ without adiabatic heating. If 
the probability for an electron to escape from the shock 
transition region is given by a Dirac delta function, we 
have a bi-Maxwellian distributions which agrees with the 
simulations. 

The microphysical mechanism by which the collision- 
less shocks are sustained over the course of the simulation 
long after the initial shock formation stage is an impor- 
tant issue. We have found that this is naturally explained 
by a modified two-stream instability due to the incom- 
ing and refiecting ions in the shock transition region-the 
unstable mode being therefore along the k^ axis (shock 
normal direction) as seen from the spectral analysis of Ex 



field. The maximum growth rate, 7 = 4 x 10 



"pej 



fror 



the modified two-stream instability is enough to excite 
the observed electric field fluctuation during the shock 
transit time(~ TiQ/tOpe < I/7). We also found that the 
fluctuating fleld is responsible for the observed entropy 
generation throughout the downstream region 
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A. Appendix - Moving Wall 

Particles will interact with the wall moving in the posi- 
tive X direction at some time between t„ and tn+i and not 
exactly at t^+i- We trace back both the particles position 
in space and the walls position in space, for each particle, 
and then calculate the correct reflection velocity used to 
determine the particles final and current locations. The 
meeting location between wall and particle can be found 
by determining the time increment Aij, = tn+i — tmeet 
at which the particle collides with the wall using 



VplWx — X'lijdll V'lijdIlLWx 



(23) 



where Xp is the particle location in x at t„+i, Vp is the 
particle's x velocity at i„+i, and x^aii is given by L^ + 
v^aii^t. Once Atj, is determined the particle's position is 
advanced backward by an increment VpAt^, the current 
deposited is removed and the proper reflection velocity is 
determined. 

We determine the rebound velocities of the particles 
by Lorentz boosting into the frame co-moving with the 



wall. In this frame, the x component of the velocity is 
simply reversed and the y and z components of the ve- 
locity are unchanged after a collision. We then transform 
back to the simulation frame and determine the new ve- 
locities by equating the total momentum 4 vectors in each 
frame. The particles are then advanced by v{At — At^) 
and their momenta are calculated using the new Lorentz 
factor calculated from their rebound velocities. 
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